Priority Effects
K. Carter Wynne
# PROLOG ##############################################################################
# PROJECT: Priority Effects Experiment
# PURPOSE: Determine how species arrival via seed additions affects the community assembly of restored prairie plant communities
#
# AUTHOR: Katherine Wynne (Wynnekat@msu.edu)
# COLLAB: L. Sullivan
# CREATED: 06 December 2022
#
# FILES: 1) inner_cover_summer_2022.xlsx - vegetative cover taken in the Summer (end of June)
# 2) inner_cover_fall _2022.xlsx - vegetative cover taken at peak biomass (end ofAugust)
# 3) Checklist of Missouri Flora.xlsx - contains the information about all plants found in MO
# 4) Cleaned_PE_Detailed_Species_List_2022.xlsx - dataset containing the information about only plants found in the study
#
#LAST UPDATED: 06 September 2023
# PROLOG ##############################################################################Treatment information & Setup
NON = None (no species seeded)
SIM = Simulataneous (36 species seeded on March 22nd, 2021)
LE = Lumped early (18 early-dispersing species seeded on March 23rd, 2021 then 18 late-dispersing species seeded on September 5th, 2021)
LL = Lumped late (18 late-dispersing species seeded on March 23rd, 2021 then 18 early-dispersing species seeded on September 5th, 2021)
NAT = Natural; species added according to peak dispersal date
- May 24th 2021 (VIOPED, PACPLA)
- June 4th 2021 (SISCAM, SPHOBT, CORLAN)
- June 22nd 2021 (LOBSPI, TRAOHI, CARBUS, HEURIC)
- July 11th 2021 (CORPAL, DODMEA)
- July 25th 2021 (KOEMAC, AMOCAN)
- August 2nd 2021 (LINSUL, MELVIR)
- August 21st 2021 (DALCAN, DALPUR, ACHMIL)
- September 17th 2021 (CROSAG)
- October 4th 2021 (CHAFAS, MONFIS, RATPIN, SPOHET, RUDHIR)
- October 18th 2021 (BIDARI, HELMOL, LESCAP)
- November 1st 2021 (PENDIG, ERYYUC, HYPPUN, SORNUT, SCHSCO, LIAPYC)
- November 19th 2021 (CORTRI, PYCTEN, SOLRIG)
Load libraries
Import Datasets
Dataset cleaning
2021
### Remove unnecessary columns
#### Get rid of the notes, unknown, and senensced columns
### Fall
inner_cover_fall_2021 <- inner_cover_fall_2021[,-c(6,7,8)]
#### change log
# - senesced PLASPP changed to PLAVIR since PLALAN and PLAMAJ do not senesce until Oct. (Sep 6, 2023)
# - JUNSPP changed to JUNINT since JUNINT has been observed in 2022 and 2023 (Sep 6, 2023)
### Get rid of EMP
### Fall
inner_cover_fall_2021 <- subset(inner_cover_fall_2021, Treatment != "EMP")2022
2023
# Make a unique identifier for year and season
### 2021
#### Season
inner_cover_fall_2021$Season <- rep("Fall", nrow=(inner_cover_fall_2021))
#### Year
inner_cover_fall_2021$Year <- rep("2021", nrow=(inner_cover_fall_2021))
### 2022
# Already has the identifiers
inner_cover_fall_2022$Year <- as.factor(inner_cover_fall_2022$Year)
inner_cover_summer_2022$Year <- as.factor(inner_cover_summer_2022$Year)
### 2023
#### Season
inner_cover_summer_2023$Season <- rep("Summer", nrow=(inner_cover_summer_2023))
inner_cover_fall_2023$Season <- rep("Fall", nrow=(inner_cover_fall_2023))
#### Year
inner_cover_fall_2023$Year <- rep("2023", nrow=(inner_cover_fall_2023))
inner_cover_fall_2023$Year <- as.factor(inner_cover_fall_2023$Year)
inner_cover_summer_2023$Year <- rep("2023", nrow=(inner_cover_summer_2023))
inner_cover_summer_2023$Year <- as.factor(inner_cover_summer_2023$Year)# Join years together
full_2021_data <- inner_cover_fall_2021
full_2022_data <- full_join(inner_cover_summer_2022,inner_cover_fall_2022)## Joining with `by = join_by(Season, Year, Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name)`full_2023_data <- full_join(inner_cover_summer_2023,inner_cover_fall_2023)## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`# Join all years together
### First 2021 and 2022
full_21_22_data <- full_join(full_2021_data, full_2022_data)## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`### Then 2021 and 2022 with 2023
full_year_data <- full_join(full_21_22_data, full_2023_data)## Joining with `by = join_by(Block, Treatment, SPP6, Percent_Cover,
## Scientific_Name, Season, Year)`# lump certain taxa together
### Lump the Carex for now (*****Ask Lauren about what to do about when the CARSPP > CARFRA or CARBRE?***)
## MELOFF + MELALB -> MELSPP
## LEPVIR -> LEPSPP
## CARFRA -> CARSPP
## CARBRE -> CARSPP
for(i in 1:nrow(full_year_data)) {
if(full_year_data[i,3] == "CARFRA"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "CARBRE"){full_year_data [i,3] <- "CARSPP"}
if(full_year_data[i,3] == "MELOFF"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "MELALB"){full_year_data [i,3] <- "MELSPP"}
if(full_year_data[i,3] == "LEPSPP"){full_year_data [i,3] <- "LEPVIR"}
if(full_year_data[i,3] == "ECHSPP"){full_year_data [i,3] <- "ECHCRU"}
# Make all cover data that was less than 1 = to 1 (to make the data discrete)
if(full_year_data[i,4] < 1) {full_year_data [i,4] <- 1}
}inner_cover_max <- full_year_data %>%
group_by(Block, Treatment, Year, SPP6) %>%
summarise(max_cover=max(Percent_Cover))## `summarise()` has grouped output by 'Block', 'Treatment', 'Year'. You can
## override using the `.groups` argument.##### Note I manually checked species that were found in both the datasets (ACHMIL, TRIREP, CORLAN) and confirmed that the code above took the highest cover value for a species seen twice
### Remove Bare
inner_cover_max_only <- subset(inner_cover_max, SPP6 != "BARE")
### Remove Litter
inner_cover_max_only <- subset(inner_cover_max_only, SPP6 != "LITTER")
head(inner_cover_max_only) ## # A tibble: 6 × 5
## # Groups: Block, Treatment, Year [1]
## Block Treatment Year SPP6 max_cover
## <dbl> <chr> <chr> <chr> <dbl>
## 1 1 LE 2021 ACAVIR 1
## 2 1 LE 2021 ACHMIL 1
## 3 1 LE 2021 AMATUB 1
## 4 1 LE 2021 AMBART 20
## 5 1 LE 2021 BARVUL 1
## 6 1 LE 2021 BROJAP 1Bare/Litter
### Bare dataset
inner_bare_cover <- subset(inner_cover_max, SPP6 == "BARE")
### Litter dataset
inner_litter_cover <- subset(inner_cover_max, SPP6 == "LITTER")
### Manipulate 2021 dataset to match 2022 and 2023
#### Bare Ground
bare_2021 <- subset(bare_cover_2021, Cover_Type != "Litter")
bare_2021 <- subset(bare_2021, Treatment != "EMP")
names(bare_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
bare_2021$Year <- rep("2021", nrow(bare_2021))
for(i in 1:nrow(bare_2021)) {
if(bare_2021[i,3] == "Bare"){bare_2021[i,3] <- "BARE"}
}
#### Litter
litter_2021 <- subset(bare_cover_2021 , Cover_Type != "Bare")
litter_2021 <- subset(litter_2021, Treatment != "EMP")
names(litter_2021) <- c("Block", "Treatment", "SPP6", "max_cover")
litter_2021$Year <- rep("2021", nrow(litter_2021))
for(i in 1:nrow(litter_2021)) {
if(bare_2021[i,3] == "Litter"){bare_2021[i,3] <- "LITTER"}
}
### Full_join datasets
inner_bare_cover <- full_join(inner_bare_cover, bare_2021)## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`inner_litter_cover <- full_join(inner_litter_cover, litter_2021)## Joining with `by = join_by(Block, Treatment, Year, SPP6, max_cover)`### Add species information to the dataset
#### N/A is N = native, A = non-native, G = Genus
inner_cover_max_only <- full_join(inner_cover_max_only, species_list_PE)## Joining with `by = join_by(SPP6)`#Remove unnecessary columns
inner_cover_max_reduced <- inner_cover_max_only[, -c(6,7,10)]Exploratory Analysis of Diversity
Summary Table of Diversity Indices
Below is a summary table showing the mean and SD for mean C value (C), species richness (SR), and Shannon diversity index (SDI) for each seeding treatment.
Analyzing diversity
Mean C
Mean C - Model assumptions
Distribution of C values for all treatments were right skewed, indicating the majority of native plants in all treatments were of low conservatism value. SIM had the most “even” distribution of C values and almost every value was represented in this treatment.
## Warning in geom_density(bins = 30, alpha = 0.4): Ignoring unknown parameters:
## `bins`
# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT
# only one extreme value in LE 2021 (likely because of DALCAN and VIOPED)
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Mean_C)## # A tibble: 5 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <dbl> <int> <dbl> <dbl> <lgl> <lgl>
## 1 NAT 2023 4 18 2.44 2.69 TRUE FALSE
## 2 NAT 2023 5 12 1.73 0.5 TRUE FALSE
## 3 NON 2023 6 19 1.99 1.18 TRUE FALSE
## 4 SIM 2022 5 20 1.83 3.13 TRUE FALSE
## 5 SIM 2023 2 27 2.41 3.42 TRUE FALSE# Mean C fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Mean_C)## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Mean_C 0.926 0.548
## 2 LE 2023 Mean_C 0.953 0.767
## 3 LL 2022 Mean_C 0.890 0.317
## 4 LL 2023 Mean_C 0.861 0.193
## 5 NAT 2022 Mean_C 0.901 0.379
## 6 NAT 2023 Mean_C 0.929 0.572
## 7 NON 2022 Mean_C 0.928 0.566
## 8 NON 2023 Mean_C 0.882 0.280
## 9 SIM 2022 Mean_C 0.816 0.0809
## 10 SIM 2023 Mean_C 0.879 0.263ggplot(data = Summary_table , aes(x= Mean_C)) +geom_histogram(bins =15) 
# Vast majority of points fall along the reference line
ggqqplot(Summary_table , "Mean_C", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
Mean C - Model fitting
Mean C values were significantly higher in all treatments compared to the NON treatment.
#Dropped the mixed model due to singular fit -> the variance of the random effect block is estimated close to zero and does not further inform the data.
Summary_table$Block <- as.factor(Summary_table$Block)
# Normal interaction model
#Mean_C.mod.lmer.interaction <- lmer(Mean_C~Treatment*Year+(1|Block), data =Summary_table)
#AIC(Mean_C.mod.lmer.interaction)
# Normal - additive model (has the lowest AIC)
Mean_C.mod.lm.additive <- lm(Mean_C~Treatment+Year, data =Summary_table)
## Summary results
summary(Mean_C.mod.lm.additive)##
## Call:
## lm(formula = Mean_C ~ Treatment + Year, data = Summary_table)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07124 -0.32383 -0.07948 0.28266 1.12107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.6602 0.1608 10.327 2.16e-14 ***
## TreatmentLL 0.5149 0.2075 2.481 0.016251 *
## TreatmentNAT -0.3549 0.2075 -1.710 0.093003 .
## TreatmentNON -1.3445 0.2075 -6.479 2.89e-08 ***
## TreatmentSIM 0.7983 0.2075 3.847 0.000318 ***
## Year2023 0.2660 0.1313 2.026 0.047691 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5084 on 54 degrees of freedom
## Multiple R-squared: 0.7134, Adjusted R-squared: 0.6869
## F-statistic: 26.88 on 5 and 54 DF, p-value: 1.544e-13## AIC test
AIC(Mean_C.mod.lm.additive)## [1] 96.76091## ANOVA type 3
anova(Mean_C.mod.lm.additive)## Analysis of Variance Table
##
## Response: Mean_C
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 4 33.674 8.4186 32.5772 8.088e-14 ***
## Year 1 1.061 1.0610 4.1057 0.04769 *
## Residuals 54 13.955 0.2584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## I removed the random effect block because it had no effect on the model fit and explained no additional variance. (chose the most parsimonius model)
# r.squaredGLMM(Mean_C.mod.lmer.additive)
## Inclusion of the random effect marginally increased model fit
est.lm.treatment <- emmeans::emmeans(Mean_C.mod.lm.additive, "Treatment")
#Posthoc test
## Treatment
pairs(est.lm.treatment, adjust = "tukey")## contrast estimate SE df t.ratio p.value
## LE - LL -0.515 0.208 54 -2.481 0.1102
## LE - NAT 0.355 0.208 54 1.710 0.4366
## LE - NON 1.345 0.208 54 6.479 <.0001
## LE - SIM -0.798 0.208 54 -3.847 0.0028
## LL - NAT 0.870 0.208 54 4.191 0.0009
## LL - NON 1.859 0.208 54 8.960 <.0001
## LL - SIM -0.283 0.208 54 -1.366 0.6519
## NAT - NON 0.990 0.208 54 4.769 0.0001
## NAT - SIM -1.153 0.208 54 -5.557 <.0001
## NON - SIM -2.143 0.208 54 -10.325 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 5 estimates# Residual plot looks decent
plot(Mean_C.mod.lm.additive)
Mean C - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.c <- Summary_table_mean$C_Mean - Summary_table_mean$C_SD
upper.c <- Summary_table_mean$C_Mean + Summary_table_mean$C_SD
pd <- position_dodge(width = 0.8)
C.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = C_Mean, color = Treatment))+
geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.c , ymax = upper.c), width = 0.5, position = pd) +
labs(y = "Mean C", x = "Year" )+
theme_classic() +
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))
C.plot
Species Richness
SIM treatments had the greatest species richness compared to all other treatments, while NON and NAT had the lowest.
SR - Model assumptions
# A couple of outliers in the NAT treatment but likely alright since the variablity in outcomes seems to be a characteristic of NAT
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Species_richness)## # A tibble: 5 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <fct> <int> <dbl> <dbl> <lgl> <lgl>
## 1 LE 2023 3 15 2.10 2.2 TRUE FALSE
## 2 LL 2022 2 25 2.31 2.06 TRUE FALSE
## 3 NAT 2022 4 22 2.18 2.31 TRUE TRUE
## 4 NAT 2022 6 15 2.10 1.09 TRUE FALSE
## 5 NON 2023 6 19 1.99 1.18 TRUE FALSE# Species richness fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Species_richness)## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Species_richness 0.996 0.998
## 2 LE 2023 Species_richness 0.847 0.149
## 3 LL 2022 Species_richness 0.814 0.0778
## 4 LL 2023 Species_richness 0.825 0.0969
## 5 NAT 2022 Species_richness 0.836 0.122
## 6 NAT 2023 Species_richness 0.884 0.286
## 7 NON 2022 Species_richness 0.960 0.817
## 8 NON 2023 Species_richness 0.955 0.784
## 9 SIM 2022 Species_richness 0.893 0.332
## 10 SIM 2023 Species_richness 0.873 0.238ggplot(data = Summary_table, aes(x= Species_richness)) +geom_histogram(bins =15) 
# Vast majority of points fall along the reference line
ggqqplot(Summary_table, "Species_richness", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
SR - Model fitting
# Model Fitting
## Despite being count data, the distribution of Species richness is remarkabley balanced and not skewed. Model comparison using AIC reveals that the Normal model works better than the Poisson model. Model results also differ likely because of under dispersion. The more complicated Quasipoisson model that accounts for under dispersion produces similar results to the Normal mixed model. Therefore, I'm going to go with the less complicated Normal model.
SR.mod.normal <- lmer(Species_richness~Treatment+Year+(1|Block), data =Summary_table)
summary(SR.mod.normal)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Species_richness ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 283.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3018 -0.7002 0.0875 0.6078 1.6946
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 2.378 1.542
## Residual 7.355 2.712
## Number of obs: 60, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 19.9333 1.0639 20.8212 18.736 1.64e-14 ***
## TreatmentLL 0.3333 1.1072 49.0000 0.301 0.764634
## TreatmentNAT -1.7500 1.1072 49.0000 -1.581 0.120398
## TreatmentNON -4.3333 1.1072 49.0000 -3.914 0.000280 ***
## TreatmentSIM 4.2500 1.1072 49.0000 3.839 0.000355 ***
## Year2023 -1.7000 0.7002 49.0000 -2.428 0.018911 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.520
## TreatmntNAT -0.520 0.500
## TreatmntNON -0.520 0.500 0.500
## TreatmntSIM -0.520 0.500 0.500 0.500
## Year2023 -0.329 0.000 0.000 0.000 0.000## Treatment and Year significant
## Interaction between Treatment and Year not significant
anova(SR.mod.normal, type = 3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 474.77 118.69 4 49 16.1381 1.704e-08 ***
## Year 43.35 43.35 1 49 5.8941 0.01891 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Marginal R^2 (variance explained by ONLY fixed effects) = 0.4547589
## Conditional R^2 (variance explained by the entire model) = 0.5554578
## Inclusion of the random effect increased model fit
r.squaredGLMM(SR.mod.normal)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.## R2m R2c
## [1,] 0.4743091 0.6027586#Posthoc test
est.lmer.SR.treatment <- emmeans::emmeans(SR.mod.normal, "Treatment")
## Treatment
pairs(est.lmer.SR.treatment , adjust = "tukey")## contrast estimate SE df t.ratio p.value
## LE - LL -0.333 1.11 49 -0.301 0.9981
## LE - NAT 1.750 1.11 49 1.581 0.5166
## LE - NON 4.333 1.11 49 3.914 0.0025
## LE - SIM -4.250 1.11 49 -3.839 0.0031
## LL - NAT 2.083 1.11 49 1.882 0.3406
## LL - NON 4.667 1.11 49 4.215 0.0010
## LL - SIM -3.917 1.11 49 -3.538 0.0076
## NAT - NON 2.583 1.11 49 2.333 0.1519
## NAT - SIM -6.000 1.11 49 -5.419 <.0001
## NON - SIM -8.583 1.11 49 -7.753 <.0001
##
## Results are averaged over the levels of: Year
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates# Poisson Model
# SR.mod.pois <- glmer(Species_richness~Treatment+Year+(1|Block),family = "poisson", data =Summary_table)
# summary(SR.mod.pois)
# Anova(SR.mod.pois, type = 3)
# Check dispersion
## Underdispersed! 0.659
## Interaction between treatment and year are not significant
# Dispersion_glmer(SR.mod.pois)
# Not accounting for the underdispersion mattered a lot! .
## Used a quassipoisson model instead and year is significant
#m6 <-glmmPQL(Species_richness~Treatment+Year,
# random = ~ 1 | Block,
# family = quasipoisson(link='log'),
# data = Summary_table)
# summary(m6)
# Anova(m6, type = 3)# Residual plot looks decent except for one large outlier
plot(SR.mod.normal)
SR - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.SR.c <- Summary_table_mean$SR_Mean - Summary_table_mean$SR_SD
upper.SR.c <- Summary_table_mean$SR_Mean + Summary_table_mean$SR_SD
pd <- position_dodge(width = 0.8)
SR.plot1 <- ggplot(data = Summary_table_mean, aes(x = Year, y = SR_Mean, color = Treatment))+
geom_point(data = Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.SR.c , ymax = upper.SR.c), width = 0.5, position = pd) +
labs(y = "Total Species Richness", x = "Year" )+
theme_classic() +
ylim(10,30)+
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
theme(legend.position="bottom")
SR.plot <- SR.plot1 +
theme(legend.position="none")
leg <- get_legend(SR.plot1)Shannon Diversity Index
SDI - Model assumptions
SIM and LL treatments had the greatest Shannon diversity index value while NON had the lowest. NAT and LE appear comparable.
# Two extreme outliers (one in NAT 2021 and another in SIM 2021)
Summary_table %>%
group_by(Treatment, Year) %>%
identify_outliers(Shannon)## # A tibble: 1 × 8
## Treatment Year Block Species_richness Shannon Mean_C is.outlier is.extreme
## <chr> <chr> <fct> <int> <dbl> <dbl> <lgl> <lgl>
## 1 NON 2022 3 11 1.55 0.429 TRUE FALSE# Except for one point, Shannon diversity fits a normal distribution
Summary_table %>%
group_by(Treatment, Year) %>%
shapiro_test(Shannon)## # A tibble: 10 × 5
## Treatment Year variable statistic p
## <chr> <chr> <chr> <dbl> <dbl>
## 1 LE 2022 Shannon 0.971 0.901
## 2 LE 2023 Shannon 0.967 0.873
## 3 LL 2022 Shannon 0.977 0.936
## 4 LL 2023 Shannon 0.854 0.169
## 5 NAT 2022 Shannon 0.950 0.743
## 6 NAT 2023 Shannon 0.931 0.586
## 7 NON 2022 Shannon 0.848 0.151
## 8 NON 2023 Shannon 0.950 0.739
## 9 SIM 2022 Shannon 0.872 0.235
## 10 SIM 2023 Shannon 0.940 0.663ggplot(data = Summary_table, aes(x= Shannon)) +geom_histogram(bins =15) 
# Vast majority of points fall along the reference line
ggqqplot(Summary_table, "Shannon", ggtheme = theme_bw()) +
facet_grid(Year~ Treatment, labeller = "label_both")
SDI - Model fitting
Shannon.lmer.mod <- lmer(Shannon~Treatment+Year+(1|Block), data =Summary_table)
summary(Shannon.lmer.mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Shannon ~ Treatment + Year + (1 | Block)
## Data: Summary_table
##
## REML criterion at convergence: 2.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.60837 -0.47607 0.08046 0.51127 2.04855
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.007995 0.08942
## Residual 0.042202 0.20543
## Number of obs: 60, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.16123 0.07452 28.51199 29.003 < 2e-16 ***
## TreatmentLL 0.05159 0.08387 49.00000 0.615 0.54133
## TreatmentNAT -0.12219 0.08387 49.00000 -1.457 0.15152
## TreatmentNON -0.34644 0.08387 49.00000 -4.131 0.00014 ***
## TreatmentSIM 0.12917 0.08387 49.00000 1.540 0.12994
## Year2023 -0.01536 0.05304 49.00000 -0.290 0.77333
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtNON TrtSIM
## TreatmentLL -0.563
## TreatmntNAT -0.563 0.500
## TreatmntNON -0.563 0.500 0.500
## TreatmntSIM -0.563 0.500 0.500 0.500
## Year2023 -0.356 0.000 0.000 0.000 0.000anova(Shannon.lmer.mod, type = 3)## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment 1.65266 0.41316 4 49 9.7902 6.691e-06 ***
## Year 0.00354 0.00354 1 49 0.0839 0.7733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# No significant interaction between treatment and year
## Marginal R^2 (variance explained by ONLY fixed effects) = 0.330704
## Conditional R^2 (variance explained by the entire model) = 0.3844412
## Inclusion of the random effect increased model fit
r.squaredGLMM(Shannon.lmer.mod)## R2m R2c
## [1,] 0.358652 0.4608049# Post-Hoc comparison
pairs(emmeans::emmeans(Shannon.lmer.mod, "Treatment"))## contrast estimate SE df t.ratio p.value
## LE - LL -0.0516 0.0839 49 -0.615 0.9720
## LE - NAT 0.1222 0.0839 49 1.457 0.5947
## LE - NON 0.3464 0.0839 49 4.131 0.0013
## LE - SIM -0.1292 0.0839 49 -1.540 0.5420
## LL - NAT 0.1738 0.0839 49 2.072 0.2486
## LL - NON 0.3980 0.0839 49 4.746 0.0002
## LL - SIM -0.0776 0.0839 49 -0.925 0.8858
## NAT - NON 0.2242 0.0839 49 2.674 0.0727
## NAT - SIM -0.2514 0.0839 49 -2.997 0.0331
## NON - SIM -0.4756 0.0839 49 -5.671 <.0001
##
## Results are averaged over the levels of: Year
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 5 estimates# A little heteroscedastic (bowtie shaped but probably ok)
plot(Shannon.lmer.mod)
SDI - Figures
Summary_table_mean$Treatment <- factor(Summary_table_mean$Treatment, levels=c("NON", "NAT", "LE", "LL", "SIM"))
lower.SDI.c <- Summary_table_mean$SDI_Mean - Summary_table_mean$SDI_SD
upper.SDI.c <- Summary_table_mean$SDI_Mean + Summary_table_mean$SDI_SD
pd <- position_dodge(width = 0.8)
SDI.plot <- ggplot(data = Summary_table_mean, aes(x = Year, y = SDI_Mean, color = Treatment))+
geom_point(data =Summary_table_mean, size = 3, aes(shape = Treatment, fill = Treatment), position = pd)+
geom_errorbar(aes(ymin = lower.SDI.c , ymax = upper.SDI.c), width = 0.5, position = pd) +
labs(y = "SDI", x = "Year" )+
theme_classic() +
theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), legend.position="none")+
scale_color_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_fill_manual(name = "Treatment", values=c("#E69F00", "#F0E442", "#009E73", "#56B4E9", "#0072B2"), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))+
scale_shape_manual(name = "Treatment", values=c(21, 22,23, 24, 25), breaks=c("NON", "NAT", "LE", "LL", "SIM"), labels = c("NON", "NAT", "LE", "LL", "SIM"))
SDI.plot 
Plot Panel
Comparing seeded cover (in the inner 1 m^2 sampling area) for each treatment
note some of the species used in the seed mix were present in the seed bank at my study site (e.g., RUDHIR and PENDIG). Additionally, some spillover occurred with seeded species seeding into other treatments like NON (e.g., BIDARI)
Total Seeded Cover
TSC - Model fitting
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment * Year,
## family = binomial, data = Total_Cover_Tot)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -6.4760 -2.0224 -0.4559 1.5998 4.5960
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.11998 0.08881 -23.870 < 2e-16 ***
## TreatmentLL 0.94324 0.10412 9.059 < 2e-16 ***
## TreatmentNAT -0.38112 0.13799 -2.762 0.00575 **
## TreatmentSIM 1.06270 0.10397 10.221 < 2e-16 ***
## Year2023 0.62783 0.11282 5.565 2.63e-08 ***
## TreatmentLL:Year2023 -0.41302 0.13685 -3.018 0.00254 **
## TreatmentNAT:Year2023 -0.36533 0.18174 -2.010 0.04441 *
## TreatmentSIM:Year2023 -0.37777 0.13662 -2.765 0.00569 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 837.79 on 47 degrees of freedom
## Residual deviance: 293.23 on 40 degrees of freedom
## AIC: 557.58
##
## Number of Fisher Scoring iterations: 5## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 272.445 3 < 2.2e-16 ***
## Year 32.035 1 1.514e-08 ***
## Treatment:Year 10.305 3 0.01614 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE Year2022 / LL Year2022 0.389 0.0405 Inf 1 -9.059 <.0001
## LE Year2022 / NAT Year2022 1.464 0.2020 Inf 1 2.762 0.1048
## LE Year2022 / SIM Year2022 0.346 0.0359 Inf 1 -10.221 <.0001
## LE Year2022 / LE Year2023 0.534 0.0602 Inf 1 -5.565 <.0001
## LE Year2022 / LL Year2023 0.314 0.0328 Inf 1 -11.076 <.0001
## LE Year2022 / NAT Year2023 1.126 0.1469 Inf 1 0.909 0.9853
## LE Year2022 / SIM Year2023 0.269 0.0281 Inf 1 -12.574 <.0001
## LL Year2022 / NAT Year2022 3.760 0.4466 Inf 1 11.150 <.0001
## LL Year2022 / SIM Year2022 0.887 0.0680 Inf 1 -1.558 0.7751
## LL Year2022 / LE Year2023 1.371 0.1210 Inf 1 3.572 0.0085
## LL Year2022 / LL Year2023 0.807 0.0625 Inf 1 -2.774 0.1016
## LL Year2022 / NAT Year2023 2.892 0.3181 Inf 1 9.654 <.0001
## LL Year2022 / SIM Year2023 0.691 0.0534 Inf 1 -4.784 <.0001
## NAT Year2022 / SIM Year2022 0.236 0.0280 Inf 1 -12.169 <.0001
## NAT Year2022 / LE Year2023 0.365 0.0461 Inf 1 -7.977 <.0001
## NAT Year2022 / LL Year2023 0.215 0.0256 Inf 1 -12.917 <.0001
## NAT Year2022 / NAT Year2023 0.769 0.1096 Inf 1 -1.842 0.5909
## NAT Year2022 / SIM Year2023 0.184 0.0219 Inf 1 -14.231 <.0001
## SIM Year2022 / LE Year2023 1.545 0.1361 Inf 1 4.935 <.0001
## SIM Year2022 / LL Year2023 0.909 0.0702 Inf 1 -1.234 0.9218
## SIM Year2022 / NAT Year2023 3.259 0.3580 Inf 1 10.753 <.0001
## SIM Year2022 / SIM Year2023 0.779 0.0600 Inf 1 -3.246 0.0258
## LE Year2023 / LL Year2023 0.588 0.0523 Inf 1 -5.971 <.0001
## LE Year2023 / NAT Year2023 2.109 0.2495 Inf 1 6.312 <.0001
## LE Year2023 / SIM Year2023 0.504 0.0447 Inf 1 -7.729 <.0001
## LL Year2023 / NAT Year2023 3.585 0.3958 Inf 1 11.563 <.0001
## LL Year2023 / SIM Year2023 0.857 0.0667 Inf 1 -1.988 0.4901
## NAT Year2023 / SIM Year2023 0.239 0.0264 Inf 1 -12.981 <.0001
##
## P value adjustment: tukey method for comparing a family of 8 estimates
## Tests are performed on the log odds ratio scale
TSC - Figures
Early season seeded cover
ESC - Model fitting
##
## Call:
## glm(formula = cbind(seeded_cover, tot_cover) ~ Treatment + Year,
## family = binomial, data = Total_Cover_Early)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8127 -1.5620 -0.8626 0.7913 5.1753
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9526 0.1056 -27.955 < 2e-16 ***
## TreatmentLL -3.3561 0.3411 -9.838 < 2e-16 ***
## TreatmentNAT -2.0984 0.1995 -10.519 < 2e-16 ***
## TreatmentSIM -0.3623 0.1070 -3.387 0.000707 ***
## Year2023 1.0821 0.1106 9.780 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 734.86 on 47 degrees of freedom
## Residual deviance: 247.79 on 43 degrees of freedom
## AIC: 381.05
##
## Number of Fisher Scoring iterations: 6## Analysis of Deviance Table (Type III tests)
##
## Response: cbind(seeded_cover, tot_cover)
## LR Chisq Df Pr(>Chisq)
## Treatment 373.50 3 < 2.2e-16 ***
## Year 106.69 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 28.6759 9.7820 Inf 1 9.838 <.0001
## LE / NAT 8.1535 1.6265 Inf 1 10.519 <.0001
## LE / SIM 1.4366 0.1537 Inf 1 3.387 0.0039
## LL / NAT 0.2843 0.1089 Inf 1 -3.284 0.0056
## LL / SIM 0.0501 0.0172 Inf 1 -8.706 <.0001
## NAT / SIM 0.1762 0.0360 Inf 1 -8.502 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale
ESC - Figures
Fall dispersing species only
LSC - Model fitting
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(seeded_cover, tot_cover) ~ Treatment + Year + (1 | Block)
## Data: Total_Cover_Late
##
## AIC BIC logLik deviance df.resid
## 605.3 616.5 -296.6 593.3 42
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.3147 -1.7714 -0.6057 1.2574 6.5425
##
## Random effects:
## Groups Name Variance Std.Dev.
## Block (Intercept) 0.008129 0.09016
## Number of obs: 48, groups: Block, 6
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.74863 0.09364 -29.352 < 2e-16 ***
## TreatmentLL 1.59358 0.09078 17.554 < 2e-16 ***
## TreatmentNAT 0.16919 0.11161 1.516 0.12953
## TreatmentSIM 1.55200 0.09191 16.885 < 2e-16 ***
## Year2023 0.13908 0.05056 2.751 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLL TrtNAT TrtSIM
## TreatmentLL -0.796
## TreatmntNAT -0.643 0.664
## TreatmntSIM -0.789 0.807 0.656
## Year2023 -0.279 0.018 -0.002 0.024## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: cbind(seeded_cover, tot_cover)
## Chisq Df Pr(>Chisq)
## (Intercept) 861.5598 1 < 2e-16 ***
## Treatment 566.5093 3 < 2e-16 ***
## Year 7.5683 1 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## contrast odds.ratio SE df null z.ratio p.value
## LE / LL 0.203 0.0184 Inf 1 -17.554 <.0001
## LE / NAT 0.844 0.0942 Inf 1 -1.516 0.4278
## LE / SIM 0.212 0.0195 Inf 1 -16.885 <.0001
## LL / NAT 4.155 0.3536 Inf 1 16.738 <.0001
## LL / SIM 1.042 0.0592 Inf 1 0.732 0.8841
## NAT / SIM 0.251 0.0216 Inf 1 -16.031 <.0001
##
## Results are averaged over the levels of: Year
## P value adjustment: tukey method for comparing a family of 4 estimates
## Tests are performed on the log odds ratio scale
LSC - Figures
Plot Panel
Community Analysis
RDA Ordination
Year: 2022
Year: 2023
### fit a distance-based RDA
#### Allows you to do a constrained ordination using non-Euclidean distance measures such as Bray-curtis dissimilarity (which is fitting for our data on plant communities)
#### it works by: 1) making the distance matrix based on your desired measure, 2) doing a PCoA, 3) taking the PCoA eigenvectors and putting them into an RDA.
### Use the capscale function (not the dbRDA since it doesn't provide species scores)
dbRDA <- capscale(inner.max.2023 ~ Treatment+Middle, data = env.data.2023, dist="bray")
## dbRDA results
dbRDA## Call: capscale(formula = inner.max.2023 ~ Treatment + Middle, data =
## env.data.2023, distance = "bray")
##
## Inertia Proportion Rank
## Total 8.39021 1.00000
## Constrained 2.96793 0.35374 5
## Unconstrained 5.65598 0.67412 24
## Imaginary -0.23369 -0.02785 5
## Inertia is squared Bray distance
## Species scores projected from 'inner.max.2023'
##
## Eigenvalues for constrained axes:
## CAP1 CAP2 CAP3 CAP4 CAP5
## 1.3067 0.6666 0.4553 0.3562 0.1831
##
## Eigenvalues for unconstrained axes:
## MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
## 0.7930 0.6676 0.5657 0.4643 0.4251 0.3832 0.3259 0.3169
## (Showing 8 of 24 unconstrained eigenvalues)## PERMANOVA results
anova(dbRDA, step = 1000, by = "terms", permu = 200)## Permutation test for capscale under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## Model: capscale(formula = inner.max.2023 ~ Treatment + Middle, data = env.data.2023, distance = "bray")
## Df SumOfSqs F Pr(>F)
## Treatment 4 2.4677 2.6178 0.001 ***
## Middle 1 0.5002 2.1227 0.012 *
## Residual 24 5.6560
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot2 <- ordiplot(dbRDA, choices = c(1,2))
sites.long2 <- sites.long(plot2, env.data = env.data.2023)
head(sites.long2)## Block Treatment Year Bare_cover Above Ground Middle axis1 axis2
## 1 1 LE 2023 3 1852.0 114.0 1305.5 -0.46435288 -1.4255688
## 2 1 LL 2023 3 1840.0 59.0 1363.0 -0.36081765 1.2530781
## 3 1 NAT 2023 2 1854.5 127.0 1335.0 -0.04903448 -0.3763299
## 4 1 NON 2023 4 1806.0 43.5 671.0 1.04202886 0.8878276
## 5 1 SIM 2023 13 1819.5 87.5 1486.5 -1.18095045 -0.3935962
## 6 2 LE 2023 18 1786.0 196.5 1252.0 0.19736143 -1.4124790
## labels
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6species.long2 <- species.long(plot2)
species.long2## axis1 axis2 labels
## ACHMIL -0.058008984 -0.510208190 ACHMIL
## AGRHYE -0.090320597 -0.076520219 AGRHYE
## ALLVIN -0.062351456 -0.055423685 ALLVIN
## AMBART 0.036005771 -0.127031844 AMBART
## AMBTRI 0.255047389 0.115741012 AMBTRI
## BARVUL 0.075742234 0.080996519 BARVUL
## BIDARI -0.070199953 -0.046741848 BIDARI
## BROJAP 0.087687511 0.090800725 BROJAP
## CARSPP -0.008329482 -0.025604545 CARSPP
## CERSPP -0.156534913 -0.049741117 CERSPP
## CHAFAS -0.214205742 0.070337697 CHAFAS
## CONCAN 0.254216008 0.144906241 CONCAN
## CORLAN -0.103646306 -0.238880443 CORLAN
## CORTRI -0.372322690 0.038493999 CORTRI
## CYNDAC 0.105773316 0.050351793 CYNDAC
## DESSPP 0.046739204 0.009500142 DESSPP
## ECHCRU -0.083299566 -0.019166459 ECHCRU
## ERASPE -0.138689185 0.065048505 ERASPE
## ERIANN -0.029273463 0.071860807 ERIANN
## ERYYUC -0.318825563 0.109071207 ERYYUC
## FESARU -0.041648471 -0.062619061 FESARU
## GALAPA 0.044397969 -0.108421281 GALAPA
## GALPED -0.099754338 -0.029167951 GALPED
## GERCAR 0.124078409 0.132483643 GERCAR
## GEUCAN 0.015747167 -0.009337348 GEUCAN
## HELAUT 0.056640393 0.025594942 HELAUT
## HELMOL -0.353685656 0.164837885 HELMOL
## HORPUS 0.195494739 0.220654291 HORPUS
## IPOHED 0.105773316 0.050351793 IPOHED
## IPOLAC -0.014429565 0.060524586 IPOLAC
## IVAANN -0.154551424 -0.046742336 IVAANN
## JUNINT -0.108982302 0.070549160 JUNINT
## KOEMAC -0.008746889 -0.275074186 KOEMAC
## LEPVIR 0.167659097 -0.206725676 LEPVIR
## LESCAP -0.216729478 0.145955400 LESCAP
## LIAPYC -0.121062510 -0.042119423 LIAPYC
## LONMAC -0.127721583 -0.181574047 LONMAC
## LOTCOR -0.179556651 -0.203237739 LOTCOR
## MELSPP 0.171183400 0.083139641 MELSPP
## MONFIS -0.084054969 -0.191315664 MONFIS
## MYOVER 0.306776454 -0.235079781 MYOVER
## OXADIL 0.177099525 -0.025248593 OXADIL
## PANDIC 0.018092319 -0.046244323 PANDIC
## PENDIG 0.157182736 -0.064837910 PENDIG
## PHAARU -0.016133670 -0.026788672 PHAARU
## PLALAN -0.016041685 -0.003680926 PLALAN
## PLAMAJ 0.105361035 0.004965493 PLAMAJ
## PLAVIR 0.176181476 -0.042260049 PLAVIR
## POAPRA -0.155178506 -0.044980390 POAPRA
## PYCTEN 0.033802919 0.001637247 PYCTEN
## RANABO -0.057204692 0.004765482 RANABO
## RATPIN -0.519755867 0.358576742 RATPIN
## RUDHIR -0.158255667 0.026219980 RUDHIR
## RUMCRI 0.044999745 -0.073745584 RUMCRI
## SETFAB 0.319110355 0.170382494 SETFAB
## SETPUM 0.565939547 0.007443133 SETPUM
## SIDSPI 0.102095417 0.041813388 SIDSPI
## SILANT 0.106599037 0.050853680 SILANT
## SOLCAN 0.083506391 0.036817574 SOLCAN
## SOLCAR 0.154563253 0.069255882 SOLCAR
## SORNUT -0.458619292 0.200809834 SORNUT
## SPHOBT -0.127212271 -0.014392910 SPHOBT
## SPOHET -0.121062510 -0.042119423 SPOHET
## SYMPIL 0.078367533 -0.137094939 SYMPIL
## TAROFF 0.119137849 -0.062499146 TAROFF
## TRIPRA -0.121062510 -0.042119423 TRIPRA
## TRIREP 0.176311277 -0.181646697 TRIREP
## VERARV 0.412681869 0.226543389 VERARV
## VERURT -0.108982302 0.070549160 VERURT
## VIOPED 0.044027264 -0.176350164 VIOPEDaxis.long2 <- axis.long(dbRDA, choices = c(1,2))
axis.long2 ## axis ggplot label
## 1 1 xlab.label CAP1 (15.6%)
## 2 2 ylab.label CAP2 (7.9%)BioR.theme <- theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line("gray25"),
text = element_text(size = 12, family="Arial"),
axis.text = element_text(size = 10, colour = "gray25"),
axis.title = element_text(size = 14, colour = "gray25"),
legend.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.key = element_blank())plotgg2 <- ggplot() +
geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +
xlab(axis.long2[1, "label"]) +
ylab(axis.long2[2, "label"]) +
scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
geom_point(data=sites.long2,
aes(x=axis1, y=axis2, colour=Treatment, shape=Treatment),
size=5) +
geom_point(data=species.long2,
aes(x=axis1, y=axis2)) +
BioR.theme +
ggsci::scale_colour_npg() +
coord_fixed(ratio=1)spec.envfit <- envfit(plot2, env=inner.max.2023)
spec.data.envfit <- data.frame(r=spec.envfit$vectors$r, p=spec.envfit$vectors$pvals)
species.long2 <- species.long(plot2, spec.data=spec.data.envfit)
species.long2## r p axis1 axis2 labels
## ACHMIL 0.5072454500 0.001 -0.058008984 -0.510208190 ACHMIL
## AGRHYE 0.0639828837 0.407 -0.090320597 -0.076520219 AGRHYE
## ALLVIN 0.0377623714 0.611 -0.062351456 -0.055423685 ALLVIN
## AMBART 0.0599078832 0.441 0.036005771 -0.127031844 AMBART
## AMBTRI 0.2392522931 0.028 0.255047389 0.115741012 AMBTRI
## BARVUL 0.1303343327 0.143 0.075742234 0.080996519 BARVUL
## BIDARI 0.0168963769 0.820 -0.070199953 -0.046741848 BIDARI
## BROJAP 0.0651875300 0.398 0.087687511 0.090800725 BROJAP
## CARSPP 0.0248129138 0.744 -0.008329482 -0.025604545 CARSPP
## CERSPP 0.0287982535 0.672 -0.156534913 -0.049741117 CERSPP
## CHAFAS 0.2000449585 0.033 -0.214205742 0.070337697 CHAFAS
## CONCAN 0.2824451407 0.010 0.254216008 0.144906241 CONCAN
## CORLAN 0.2291601219 0.022 -0.103646306 -0.238880443 CORLAN
## CORTRI 0.3719211480 0.001 -0.372322690 0.038493999 CORTRI
## CYNDAC 0.1014468817 0.309 0.105773316 0.050351793 CYNDAC
## DESSPP 0.0150879971 0.912 0.046739204 0.009500142 DESSPP
## ECHCRU 0.0465040615 0.694 -0.083299566 -0.019166459 ECHCRU
## ERASPE 0.1031003723 0.251 -0.138689185 0.065048505 ERASPE
## ERIANN 0.0522202337 0.531 -0.029273463 0.071860807 ERIANN
## ERYYUC 0.2378151015 0.023 -0.318825563 0.109071207 ERYYUC
## FESARU 0.0008621138 0.991 -0.041648471 -0.062619061 FESARU
## GALAPA 0.0498040088 0.692 0.044397969 -0.108421281 GALAPA
## GALPED 0.0297373414 0.647 -0.099754338 -0.029167951 GALPED
## GERCAR 0.1152105947 0.187 0.124078409 0.132483643 GERCAR
## GEUCAN 0.0172861958 0.853 0.015747167 -0.009337348 GEUCAN
## HELAUT 0.0128013368 0.917 0.056640393 0.025594942 HELAUT
## HELMOL 0.4037999010 0.002 -0.353685656 0.164837885 HELMOL
## HORPUS 0.2497223863 0.011 0.195494739 0.220654291 HORPUS
## IPOHED 0.1014468817 0.309 0.105773316 0.050351793 IPOHED
## IPOLAC 0.0269072643 0.731 -0.014429565 0.060524586 IPOLAC
## IVAANN 0.1299775945 0.119 -0.154551424 -0.046742336 IVAANN
## JUNINT 0.0994951469 0.326 -0.108982302 0.070549160 JUNINT
## KOEMAC 0.2098574477 0.025 -0.008746889 -0.275074186 KOEMAC
## LEPVIR 0.1922656027 0.075 0.167659097 -0.206725676 LEPVIR
## LESCAP 0.2672247308 0.003 -0.216729478 0.145955400 LESCAP
## LIAPYC 0.1084649946 0.240 -0.121062510 -0.042119423 LIAPYC
## LONMAC 0.2098758530 0.038 -0.127721583 -0.181574047 LONMAC
## LOTCOR 0.2377141991 0.017 -0.179556651 -0.203237739 LOTCOR
## MELSPP 0.1105157989 0.198 0.171183400 0.083139641 MELSPP
## MONFIS 0.1949009232 0.044 -0.084054969 -0.191315664 MONFIS
## MYOVER 0.2442572475 0.024 0.306776454 -0.235079781 MYOVER
## OXADIL 0.1444608064 0.091 0.177099525 -0.025248593 OXADIL
## PANDIC 0.0021488040 0.981 0.018092319 -0.046244323 PANDIC
## PENDIG 0.0933758114 0.281 0.157182736 -0.064837910 PENDIG
## PHAARU 0.0089886654 0.966 -0.016133670 -0.026788672 PHAARU
## PLALAN 0.0001244903 0.999 -0.016041685 -0.003680926 PLALAN
## PLAMAJ 0.0481384833 0.624 0.105361035 0.004965493 PLAMAJ
## PLAVIR 0.1187269285 0.158 0.176181476 -0.042260049 PLAVIR
## POAPRA 0.1689308407 0.060 -0.155178506 -0.044980390 POAPRA
## PYCTEN 0.0148643439 0.936 0.033802919 0.001637247 PYCTEN
## RANABO 0.0176453416 0.799 -0.057204692 0.004765482 RANABO
## RATPIN 0.6577890834 0.001 -0.519755867 0.358576742 RATPIN
## RUDHIR 0.0591117238 0.429 -0.158255667 0.026219980 RUDHIR
## RUMCRI 0.0039828781 0.960 0.044999745 -0.073745584 RUMCRI
## SETFAB 0.3549839356 0.001 0.319110355 0.170382494 SETFAB
## SETPUM 0.7069156513 0.001 0.565939547 0.007443133 SETPUM
## SIDSPI 0.0813242368 0.369 0.102095417 0.041813388 SIDSPI
## SILANT 0.0406318031 0.748 0.106599037 0.050853680 SILANT
## SOLCAN 0.0370632877 0.769 0.083506391 0.036817574 SOLCAN
## SOLCAR 0.1279225355 0.124 0.154563253 0.069255882 SOLCAR
## SORNUT 0.5060084303 0.001 -0.458619292 0.200809834 SORNUT
## SPHOBT 0.0259085617 0.751 -0.127212271 -0.014392910 SPHOBT
## SPOHET 0.1084649946 0.240 -0.121062510 -0.042119423 SPOHET
## SYMPIL 0.1437469470 0.123 0.078367533 -0.137094939 SYMPIL
## TAROFF 0.0595860076 0.427 0.119137849 -0.062499146 TAROFF
## TRIPRA 0.1084649946 0.240 -0.121062510 -0.042119423 TRIPRA
## TRIREP 0.2735959378 0.015 0.176311277 -0.181646697 TRIREP
## VERARV 0.3327591785 0.007 0.412681869 0.226543389 VERARV
## VERURT 0.0994951469 0.326 -0.108982302 0.070549160 VERURT
## VIOPED 0.2171418847 0.002 0.044027264 -0.176350164 VIOPEDspecies.long3 <- species.long2[species.long2$p < 0.05, ]
species.long3## r p axis1 axis2 labels
## ACHMIL 0.5072454 0.001 -0.058008984 -0.510208190 ACHMIL
## AMBTRI 0.2392523 0.028 0.255047389 0.115741012 AMBTRI
## CHAFAS 0.2000450 0.033 -0.214205742 0.070337697 CHAFAS
## CONCAN 0.2824451 0.010 0.254216008 0.144906241 CONCAN
## CORLAN 0.2291601 0.022 -0.103646306 -0.238880443 CORLAN
## CORTRI 0.3719211 0.001 -0.372322690 0.038493999 CORTRI
## ERYYUC 0.2378151 0.023 -0.318825563 0.109071207 ERYYUC
## HELMOL 0.4037999 0.002 -0.353685656 0.164837885 HELMOL
## HORPUS 0.2497224 0.011 0.195494739 0.220654291 HORPUS
## KOEMAC 0.2098574 0.025 -0.008746889 -0.275074186 KOEMAC
## LESCAP 0.2672247 0.003 -0.216729478 0.145955400 LESCAP
## LONMAC 0.2098759 0.038 -0.127721583 -0.181574047 LONMAC
## LOTCOR 0.2377142 0.017 -0.179556651 -0.203237739 LOTCOR
## MONFIS 0.1949009 0.044 -0.084054969 -0.191315664 MONFIS
## MYOVER 0.2442572 0.024 0.306776454 -0.235079781 MYOVER
## RATPIN 0.6577891 0.001 -0.519755867 0.358576742 RATPIN
## SETFAB 0.3549839 0.001 0.319110355 0.170382494 SETFAB
## SETPUM 0.7069157 0.001 0.565939547 0.007443133 SETPUM
## SORNUT 0.5060084 0.001 -0.458619292 0.200809834 SORNUT
## TRIREP 0.2735959 0.015 0.176311277 -0.181646697 TRIREP
## VERARV 0.3327592 0.007 0.412681869 0.226543389 VERARV
## VIOPED 0.2171419 0.002 0.044027264 -0.176350164 VIOPEDspecies.long3 <- species.long3 %>%
rename(SPP6 = "labels")
species.long3 <- left_join(species.long3, species_list_PE, by = "SPP6")vectors.envfit <- envfit(plot2, env=env.data.2023)
vectors.long3 <- vectorfit.long(vectors.envfit)
vectors.long3## vector axis1 axis2 r p
## Block Block 0.8032309 0.5956678 0.02300445 0.729
## Bare_cover Bare_cover 0.9653136 -0.2610932 0.01682397 0.806
## Above Above -0.4497822 -0.8931382 0.01865029 0.786
## Ground Ground -0.8836873 -0.4680778 0.06824846 0.409
## Middle Middle -0.8443056 -0.5358620 0.25754975 0.019ordiplot(dbRDA, choices=c(1,2))
Treatment.ellipses <- ordiellipse(plot2, groups=env.data.2023$Treatment, display="sites", kind="sd")
Treatment.ellipses.long2 <- ordiellipse.long(Treatment.ellipses, grouping.name="Treatment")
All
PERMANOVA
PERMANOVA analysis revealed that treatment did have a significant effect on plant community composition. Treatment also explained ~ 20% of the variation in the plant community. Post-hoc analysis shows that LE, NON, and NAT treatments were not significantly different from each other but did differ in plant community composition compared to the SIM and LL seeding treatments. Additionally, the SIM and LL treatments were only marginally different (p = 0.053). Based on the NMDS above, the LL plant community appears as a subset of the SIM plant community. Despite receiving the same pool of seeded species, species added to the plots later in the growing season did not establish as well. Barriers to establishment also appeared higher for spring/summer dispersing species compared to fall dispersing species, particularly when added without priority.
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 1000
##
## adonis2(formula = inner.max ~ Treatment + Year, data = inner_wide, permutations = 1000, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treatment 4 3.1762 0.17687 3.3089 0.000999 ***
## Year 1 1.8224 0.10148 7.5940 0.000999 ***
## Residual 54 12.9588 0.72164
## Total 59 17.9574 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## $parent_call
## [1] "inner.max ~ Treatment + Year , strata = Null , permutations 999"
##
## $LE_vs_LL
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.9839 0.14313 4.2642 0.001 ***
## Year 1 1.0450 0.15202 4.5291 0.001 ***
## Residual 21 4.8456 0.70486
## Total 23 6.8745 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.4410 0.06307 1.6582 0.039 *
## Year 1 0.9656 0.13812 3.6310 0.001 ***
## Residual 21 5.5845 0.79881
## Total 23 6.9911 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.6830 0.10287 2.8438 0.001 ***
## Year 1 0.9126 0.13746 3.8000 0.001 ***
## Residual 21 5.0433 0.75967
## Total 23 6.6389 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LE_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.7381 0.10692 2.9548 0.001 ***
## Year 1 0.9197 0.13323 3.6819 0.001 ***
## Residual 21 5.2454 0.75986
## Total 23 6.9031 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NAT
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.7043 0.10632 2.8839 0.001 ***
## Year 1 0.7917 0.11951 3.2417 0.001 ***
## Residual 21 5.1287 0.77418
## Total 23 6.6247 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 1.1298 0.17501 5.2743 0.001 ***
## Year 1 0.8276 0.12820 3.8636 0.001 ***
## Residual 21 4.4985 0.69680
## Total 23 6.4560 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $LL_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.5752 0.09413 2.6448 0.001 ***
## Year 1 0.9681 0.15842 4.4511 0.001 ***
## Residual 21 4.5672 0.74744
## Total 23 6.1105 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_NON
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.4520 0.07021 1.8176 0.026 *
## Year 1 0.7635 0.11859 3.0701 0.001 ***
## Residual 21 5.2222 0.81119
## Total 23 6.4377 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NAT_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 0.9454 0.13241 3.6145 0.001 ***
## Year 1 0.7019 0.09831 2.6836 0.003 **
## Residual 21 5.4929 0.76928
## Total 23 7.1403 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $NON_vs_SIM
## Df SumOfSqs R2 F Pr(>F)
## Treatment 1 1.2878 0.18695 5.5308 0.001 ***
## Year 1 0.7109 0.10321 3.0533 0.003 **
## Residual 21 4.8897 0.70984
## Total 23 6.8884 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"Weather data
# Precipitation
month_tot <- weather %>%
filter(Year != "2021") %>%
group_by(Month, Year) %>%
summarize(tot_month_PRP_mm = sum(Total_PRP_mm))## `summarise()` has grouped output by 'Month'. You can override using the
## `.groups` argument.year_tot <- weather %>%
filter(Year != "2021") %>%
group_by(Year) %>%
summarize(tot_PRP_mm = sum(Total_PRP_mm))
library(mosaic)## Registered S3 method overwritten by 'mosaic':
## method from
## fortify.SpatialPolygonsDataFrame ggplot2
##
## The 'mosaic' package masks several functions from core packages in order to add
## additional features. The original behavior of these functions should not be affected by this.
##
## Attaching package: 'mosaic'
##
## The following objects are masked from 'package:rstatix':
##
## cor_test, prop_test, t_test
##
## The following object is masked from 'package:lmerTest':
##
## rand
##
## The following object is masked from 'package:lme4':
##
## factorize
##
## The following object is masked from 'package:Matrix':
##
## mean
##
## The following object is masked from 'package:cowplot':
##
## theme_map
##
## The following objects are masked from 'package:goeveg':
##
## deg2rad, rad2deg
##
## The following object is masked from 'package:permute':
##
## shuffle
##
## The following objects are masked from 'package:dplyr':
##
## count, do, tally
##
## The following object is masked from 'package:purrr':
##
## cross
##
## The following object is masked from 'package:ggplot2':
##
## stat
##
## The following objects are masked from 'package:stats':
##
## binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
## quantile, sd, t.test, var
##
## The following objects are masked from 'package:base':
##
## max, mean, min, prod, range, sample, sumfavstats(year_tot$tot_PRP_mm)## min Q1 median Q3 max mean sd n missing
## 665.988 848.995 934.72 1035.304 1148.842 927.735 143.0993 10 0# Average air temperature
library(mosaic)
favstats(weather$Average_Air_Temp_C)## min Q1 median Q3 max mean sd n missing
## -21.5 4.277778 13.33333 21.55556 32.66667 12.39584 10.45545 4018 0